Hybrid lattice Boltzmann model for binary fluid mixtures 

A. Tiribocchi, 1 -!] N. Stella, 2 Q G. Gonnella, 1 ! and A. Lamura 3 '! 

1 Dipartimento di Fisica, Universita di Bari, and INFN, 
Sezione di Bari, Via Amendola 173, 70126 Bari, Italy 
2 Dipartimento di Fisica, Universita di Bari, 
Via Amendola 173, 70126 Ban, Italy 
3 Istituto Applicazioni Calcolo, CNR, 
Via Amendola 122/D, 70126 Ban, Italy 
(Dated: July 16, 2009) 

Abstract 

A hybrid lattice Boltzmann method (LBM) for binary mixtures based on the free-energy ap- 
proach is proposed. Non-ideal terms of the pressure tensor are included as a body force in the LBM 
kinetic equations, used to simulate the continuity and Navier-Stokes equations. The convection- 
diffusion equation is studied by finite difference methods. Differential operators are discretized in 
order to reduce the magnitude of spurious velocities. The algorithm has been shown to be stable 
and reproducing the correct equilibrium behavior in simple test configurations and to be Galilean 
invariant. Spurious velocities can be reduced of about an order of magnitude with respect to 
standard discretization procedure. 
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I. INTRODUCTION 



In recent years lattice Boltzmann methods (LBM) [l| have been widely used to study 
multiphase fluids fl. Examples of applications are the analysis of growth regimes in phase 
separation of binary mixtures [3| or the study of backflow effects in liquid crystal behavior 
4j. The LBM approach is well suited for dealing with complex geometries or for parallel 
implementations [1]. Moreover, in the free-energy approach 5|, the mesoscale properties of 
the fluid (interface structures, coupling with local order parameters, etc.) can be straightfor- 
wardly inserted in the LBM numerical scheme and taken under control. Due to the relevance 
of the method, it is worth to further develop LBM algorithms in order to improve numerical 
stability and accuracy, also by optimizing the use of computer resources. 

LBM dynamics is defined in terms of kinetic equations for a set of populations repre- 
senting, at each lattice site and time, the density of particles moving in one of the allowed 
directions of a given lattice. The sum over the directions i of fa is the local density of the 
fluid while the first momentum is related to the local fluid momentum. In one approach a 
forcing term is included in the kinetic equations representing the interactions between the 

. Differently, the free-energy method was originally developed 



components of the mixture 



the pressure tensor of the fluid 



tne pre 

3, y,Q 



by fixing the second moment of the populations in terms o: 
mixture 7|. It has been applied to complex fluids in Refs. 

In this paper we consider an approach similar to the one of Ref. where a free-energy 
dependent term is added as a body force in the kinetic equations. This approach traces back 



to the work of Guo et al. [l2j where a comparison with different methods to introduce the 
force is reported. With respect to the algorithm of Ref. [7| , this allows a better control of 



the continuum limit still keeping all the advantages of the free-energy method. In Ref. a 
lattice Boltzmann equation is considered for each component. Here we consider a "hybrid" 
algorithm where LBM is used to simulate Navier-Stokes equations while finite-difference 
methods are implemented to simulate the convection-diffusion equation. Such hybrid codes 
have been used for complex fluids [l^, liquid crystals [14] and thermal flows 15]. This allows 
to reduce in a relevant way the amount of required memory in systems with multi-component 
order parameters or in simulations of three-dimensional systems. 

A typical undesired effect due to discretization is the appearing of unphysical flow close to 
the interfaces. This flow, often known as spurious velocities, can severely affect the quality 
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of LBM simulations. In this work we discretize the differential operators by a procedure 
optimized for reducing the magnitude of spurious velocities, following the so-called "stencil" 
method applied in Ref. to a multiphase one-component fluid. Here we will see that this 
method allows to reduce spurious velocities of about an order of magnitude. 

The paper is organized as follows. In the next section the LBM algorithm proposed is 
described and details on the numerical implementation are given. In Section 3 results of 
simulations of test configurations are shown. We will see how spurious velocities around 
curved interfaces can be reduced applying a more general stencil to discretize derivatives. 
We will also discuss the convection of a drop under a constant force acting for a finite time 
interval. Then some conclusions will follow in Section 4. 



II. THE MODEL 



The equilibrium properties of the fluid mixture can be described by a free energy 



T = \ dr 



nTlnn + ^ 2 + ^ 4 + ^(Vy?) 2 



(1) 



where T is the temperature, n is the total density of the mixture, and ip is the scalar order 
parameter representing the concentration difference between the two components of the 
mixture. The term depending on n gives rise to the ideal gas pressure p l = nT which does 
not affect the phase behavior. The terms in (p in the free-energy density f(n, tp, T) correspond 
;o the typical expression of Ginzburg-Landau free energy used in studies of phase separation 
171 ] . The terms in the free energy can be distinguished in two parts: The polynomial terms 
describe the bulk properties of the mixture and the gradient term is related to the interfacial 
ones. 

In the bulk terms the parameter b is always positive to ensure stability while the parameter 
a can distinguish a disordered (a > 0) and an ordered (a < 0) mixture, in which the two 

components coexist with equilibrium values ±(p eq where ip eq = —g- The equilibrium 
profile between the two coexisting bulk components is 

2x 

ip(x) = V9 eg tanh(— ) (2) 

with interface width 



2k 



and surface tension 



2 2a*K 

0r= 3V— ■ (4) 



The thermodynamic functions can be obtained from the free energy by differentiation. 
The chemical potential difference between the two components is given by 

a = = clld + bip 3 — kV 2 (d. (5) 
dip 



The pressure P a p is a tensor since interfaces in the fluid can exert nonisotropic forces 19]. 
The diagonal part p Q can be obtained from Eq. ([T]) as 

p ° = + v % ~ /( "' v ' T)=p, + i^ 2 + f^ 4 _ k ^( v V) - f (v»>) 2 . (e) 

For a fluid with concentration gradients P a6 has to verify the general eqni.ibrium condition 
d a Pa/3 = [20[. A suitable choice for the pressure tensor is 

Pa(3 = PoSaf3 + Kd a lfdpif. (7) 

The hydrodynamic equations of fluids follow from the conservation laws for mass and 
momentum. For binary mixtures at constant temperature the evolution of density, velocity 
and concentration fields is described by the continuity, the Navier-Stokes and the convection- 
diffusion equations [2l|, respectively, 



d t n + d a (nu a ) = 0, (8) 

J 
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dt(nup) + d a {nu a up) = -d a P aP + d a {r](d a up + dpu a r^<9 7 u 7 ) + C<Wd 7 « 7 } 



= -dp(p l ) - (fdpfi + d a {q{d a up + dpiia 4^d 7 « 7 ) + C<W<9 7 « 7 }, (9) 

dttp + d a (ipu a ) = rvV, (io) 

where rj and ( are the shear and the bulk viscosities, T is the mobility coefficient, and d is 
the dimensionality of the system. 

Equations fl8|)- f|T0|) can be solved numerically. We use a mixed approach that consists 
of a finite difference scheme for solving Eq. (ITU]) and of a LBM approach with forcing term 
for Eqs. (jHJ) and (J9j). This has the advantage that the amount of required memory can 
be decreased so that larger systems can be simulated. In our case of study, for a two- 
dimensional model on a square lattice with nine velocities (D2Q9), this method allows to 
reduce the required memory of ~ 27%. Actually, the convection-diffusion equation could 



have also been solved on a D2Q5 lattice 22j and in this case the reduction in memory would 
have been of ~ 17%. Moreover, the spurious terms in the continuum equations found in 
previous formulations based on a free energy 0] can be avoided. 
A. The Lattice Boltzmann scheme with forcing term 



To solve Eqs. (jSJ) and (Q we use a Lattice Boltzmann scheme on a lattice of size L x x L y 

in which each site is connected to nearest and next-to-nearest neighbors. This is one of the 

simplest geometries which reproduce correctly the Navier-Stokes equations in continuum 

limit and is shown in Fig. 1. Horizontal and vertical links have length Ax and diagonal 

links \/2Ax. On each site r nine lattice velocity vectors e$ are defined. They have modulus 
Ax 

|ej| = — = c, being At^B the time step, for % = 1,2,3,4 and modulus [e» | = \/2c for 

At LB 

i = 5,6, 7, 8. Moreover, the zero velocity vector e = is defined. A set of distribution 
function {/j(r, i)} is defined on each lattice site r at each time t. 

In the LB scheme for simple fluids [1] the distribution functions evo 
step AtiB according to a single relaxation-time Boltzmann equation {23] 

At 



ve during the time 



fi(r + eiAt LB , t + At LB ) - ^(r, t) 



(11) 



where r is a relaxation parameter and /^(r, t) are the local equilibrium distribution func- 
tions. The total density n and the fluid momentum nu are defined by the following relations 



n 



(12) 



where u is the fluid velocity. The form of f^ q must be chosen so that the mass and momentum 
are locally conserved in each collision step, therefore the following relations must be satisfied 



Eur -/*)= ^E/r 

i i 



n. 



nu. 



(13) 
(14) 



Moreover, the /f 9 's need to have some symmetries so that the Navier-Stokes equations are 
reproduced in the continuum limit. A convenient choice for the local equilibrium distribution 
functions of an ideal fluid in the case of a D2Q9 model is given by a second order expansion 
in the fluid velocity u of the Mawwell-Boltzmann distribution 24| 



f- 9 (r,t) 



UOjTl 



1 + 



u 



uu :(e^ 



(15) 
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where c s = c/\f?> is the sound speed in this model, I is the unitary matrix and a suitable 
choice for the coefficients Ui is Uq = 4/9, Ui = 1/9 for i = 1 — 4, uji = 1/36 for % = 5 — 8. 
This form is such that 

fi 9 ei a e if 3 = nc 2 s 5 a/3 + nu a U(3. (16) 

i 

In order to simulate Eq. Q where a nonideal pressure tensor P a p appears, we adopt a 
LB model with a forcing term following a derivation similar to that of Ref. [Ijj]. In the case 
of Ref. 12] the model was used to study forced simple fluids while we address the case of a 
binary mixture with interaction and interface contributions. The evolution equation of the 
distribution functions becomes 



fi(r + eiAt LB , t + At lb) - /<(r, t) 



At 



LB , 



[/,(r,t)-/r(r,t)]+AW^. 



(17) 



where Fi is the forcing term to be properly determined. The equilibrium distribution func- 
tions ([T5]) are not changed except for the formal substitution u — >• u*, where u* is given 
by 



nu 



£/<• 



-FAt 



LB, 



(18) 



F being the force density acting on the fluid and u* the physical velocity. The expression 
of F for our case will be given later. The forcing term can be expressed as a power series 
at the second order in the lattice velocity 251 ] 



' B-e. t C :(e t e,-cgl) 



(19) 



where A, B, and C are functions of F. The moments of the force verify the following 
relations 



B, X)F f e i e i = ^+i[C + C a 



(20) 



and have to be consistent with the hydrodynamic equations. 

The continuum limit is obtained by using a Chapman- Enskog expansion in the Knudsen 
number e: 

/i = //° ) +6/ 4 (1) +6 a /W + ...., (21) 

d t = ed h +e 2 d t2 , (22) 



d r = ed ri , 



(23) 
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C = eCi. 



26j. The continuity and the Navier-Stokes 



F = eFi, A = eA h B = eBi 

We note that the force term is of first order in e 
equations are recovered in the following form 

dt(nu*p) + d a (nu* a u*p) = -d p (nc 2 ) + F p + d a {ri(d a Up + d p u* a )} 

in terms of the velocity u* when the following expressions for the terms A, B, C: 



(24) 



(25) 



A = 0, 



B 



At 



LB 



c 



1 



At 



LB 



(u*F + Fu*) 



(26) 



2r / V 2r 

are used. The continuum equations (jSJ) and (1251) can be also obtained by a Taylor expansion 

method. We remark that no spurious terms are present in the continuum equations except 

for a term of order u* 3 which is neglected in Eq. (1251) . Such approximation is correct as 

far as u* 2 « c 2 when the expansion (fT5l) is valid [1]. In the present formulation the 

second moment of the equilibrium distribution function (1T61) does not need to be modified 

to include the effects of the pressure tensor as in previous models based on a free energy 

[jj]. It is straightforward to show that the momentum defined in Eq. ( fTBT) corresponds to an 

average between the pre- and post-collisional values of the velocity u which is the correct 

way to calculate it when a forcing term is introduced jfj, [2^]. It is this value that appears in 

ihe continuum equations and is measured in simulations. As in the case of standard LBM 

2 

1|, the present model is characterized by the fact that £ = —r\ with shear viscosity 



d 



r\ = nc 2 s At LB 



At 



LB 



(27) 



In order to recover Eq. (Q we have to require that 



F = V(nc 2 - p 1 ) - <^V/i = -<pVfj.. 



(28) 



The last equality comes from the fact the term nc s corresponds in LBM to the ideal gas 
pressure p 1 Finally, the forcing term in Eq. ffT7|) has the form 



At 



LB 



2r 



OJi 



e, - u 



ei • u 



F 



(29) 



with u* given by Eq. ( jTBl) . 



7 



B. Numerical calculation of the forcing term 



The derivatives of the order parameter in the forcing term (1281) are calculated using a 
finite difference scheme. In particular, we have adopted a stencil representation of finite 
difference operators in the more gene ral way to ensure higher isotropy 16|, which is known 



to reduce spurious velocities 
operators are, respectively, 



27 



28( | . The schemes for the x derivative and the Laplacian 



dr)r — ~. 

Ax 



-M M 
-NON 
-MOM 



(30) 



V 2 



Ax A 



R Q R 

Q -4(Q + R) Q (31) 
R Q R 

with 2N + AM = 1 and Q + 2R = 1 to guarantee consistency between the continuous 
and discrete derivatives 16|. The subscript D in the symbols of derivatives denotes the 
discrete operator. In these schemes the central entry is referred to the lattice point where 
the derivative is computed, and the other entries are referred to the eight neighbor lattice 
sites. The discrete derivatives of the order parameter ip are computed by summing the 
values in the site and in the eight neighbors with the weights in the matrices (1301) - (1311) . The 
y derivative is computed by transposing the matrix ( |30l) . The choice of the free parameters 
N and Q is made in such a way that the spurious velocities are minimized (see next Section). 
We will refer to this the optimal choice (OC). The values N = 1/2, M = 0, Q = 1, 

and R = correspond to the standard central difference scheme denoted as SC. We will 
compare SC and OC in the following. 



C. The scheme for the convection-diffusion equation 

The convection-diffusion equation f TTTJT) is solved by using a finite difference scheme. The 

function <p(r,t) is defined on the nodes of the same lattice used for the LB scheme. The 

time is discretized in time steps AtFD with time values t n = nAtpD, n = 1,2,3,.... The 

relationship connecting the two time steps is At lb = mAtpD, being m an integer. We 
denote any discretized function at time t n on a node (x iy yj) (i = 1,2,..., L x ; j = 1, 2, L y ) 
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of the lattice by g(xi,yj,t n ) = g™. At each time step we update (p n — > ip n+1 using Eq. (flQj) 
in two successive partial steps [29]. This allows to have a better numerical stability. In the 
first step we implement the convective term using an explicit Euler algorithm 30J 



where the velocity u* comes from the solution of the LB equation. Note that the term 
d a u*™ has not been neglected since the fluid is not exactly incompressible. Indeed, the 
Navier-Stokes equation ( 1251) coming from the LBM contains some compressibility terms 
which can be anyway kept very small requiring that u* 2 « c 2 [1]. The derivatives in (|32|) 
are discretized as follows: 

d D*u*\v ~ 2Ai (33) 



duMl = ( ' + A - ^ <% < (35) 



and analogously for the y components. 

The diffusive part of Eq. ( FlOl) is implemented in the second update step using an explicit 
Euler algorithm as 

v n+l = ^n+1/2 + AtpDT [ aV 2^n+l/2 + ^2 f n _ R y 2 ( ^2^+1/2} ] (gg) 

where f n = (f n ) 3 and the operator V 2 is discretized using the form given in Eq. ( |3T|) with the 
standard choice Q — 1 and i? = 0. Other choices using a more general stencil for discretizing 
V 2 are possible though we checked that they did not provide any relevant difference. 



III. RESULTS AND DISCUSSION 



We considered several test cases in order to validate our model. We used the values 
Ax = At^B = Atp£, = 1. In the free energy we adopted the parameters —a = b = 10~ 3 , 
k = —3a corresponding to an equilibrium interface of width £ ~ 5Ax. The mobility T was 
set to 5 and the relaxation time rj At^B was 1 unless differently stated. 

We first examined the relaxation to equilibrium of a planar sharp interface on a lattice of 
size L x = L y = 64 varying r in the SC case. In all the cases the system correctly relaxes to 
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the expected profile (j2j). One example is reported in Fig. 2. In the case of a planar interface 
the fluid velocities u* decay to negligible values as it should be at equilibrium when A/z = 
and d a P a p = 0. 

We then studied a circular drop test for with interfaces not aligned with the 

lattice links. A drop with sharp interface of diameter 64Ax was placed at the center of a 
lattice of size L x = L y = 128 and let equilibrate in the SC case. Interfaces relax to the 
expected profile without deforming the drop but spurious velocities appear as it can be seen 
in the upper panel of Fig. 3 in the case with r/AtiB — 5. We then used the OC scheme to 
verify whether spurious velocities could be reduced by using a more isotropic structure for 
the discrete spatial derivatives in the forcing term (j28j) . We scanned several values of N and 
Q in order to reduce the maximum value of the velocity |«™ a J on the whole lattice. The 
optimal values are summarized in the Table I. It is interesting to note that there is a couple 
of values iV = 0.3 and Q = 2.5 which occurs more frequently. We verified that this choice is 
also effective in reducing spurious velocities even for the other values of r. For this choice of 
iV and Q the maximum velocities differ only by a small percentage from the tabled values. 
Velocities can be greatly reduced with respect to the SC case as it can be visually observed 



r/At LB 


N 


Q 


Wmax \/ c s 


0.6 


0.3 


3 


0.0001753 


0.8 


0.3 


2.5 


0.0000603 


1 


0.3 


2.5 


0.0000365 


1.2 


0.3 


2.5 


0.0000267 


5 


0.3 


2.5 


0.0000088 


10 


0.3 


2 


0.0000062 



TABLE I: Optimal values of N and Q for different values of r and the corresponding values of the 
maximum spurious velocity |it^ aa .|- 

in the lower panel of Fig. 3. 

We also tried to get an analytical estimate of the optimal values of N and Q in the 
following way. At equilibrium it holds that d a P a p = tpdpn = (cap + 3bip 3 )d/3(p — kipdpiy 2 ^) = 
0. This expression depends on the first- and third-order derivatives. By using the stencils 
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(l30]) - (l3~Tj) we get for the discretized operators the expressions 

1,. 1-2JV 



d Dx = d x + ^(Ax) 2 d 3 x + ——(Ax) 2 d x d 2 y + ... (37) 
6 2 y 



and 



so that 



V 2 D = V 2 + ^(A*) 2 (^ + + ^(Ax)^ + .., (38) 



d Dx {V 2 D ) = d x (V 2 ) + -i^xfdl - - r + ^— ^1 + — ^ (Ax) W (39) 



4 

+ 



1 1 - 2N 1-Q 



6 2 

(Ax) 2 d x di + 



■12 2 

By imposing that the error terms in the third-order derivative depending on N and Q vanish, 
we get N = 7/12 ~ 0.6 and Q = 7/6 ~ 1.2. However, this estimate does not correspond to 
the optimal results of Table I. This is due to the fact that these optimal values were found by 
considering the full dynamical problem with the whole set of equations where we minimized 
the spurious velocities. In the estimate after Eq. (TJO]) the coupling with the velocity field 
was not taken into account so that there is no a priori reason to expect the same optimal 
values for iV and Q. 

A comparison of the spurious velocities in the SC and OC cases is shown in Fig. 4. By 
using the optimal choice OC the spurious velocities can be reduced by a factor approximately 
10 with respect to the standard case SC over the whole range of r values. The stencil forms 
( |3"0]) - (!3~T|) were also applied to the model of Ref. for nonideal fluids finding a comparable 
reduction in the magnitude of spurious velocities with respect to the standard case [16| . 

We then studied the motion of an equilibrated drop of diameter 64Ax in a lattice of size 
L x = 256, L y = 128 under the effect of an external constant force that acts up to the time 
t/At LB = 500 and is then switched off. The additional force G = n(g x , 0)Ax 2 /At £B acts 
on the total density. g x is in the range [10~ 5 ,5 x 10~ 5 ] and the OC scheme is used. The 
overall system is set in motion rightwards with increasing velocity until the force G is on, 
then it moves with constant speed. The choice of g x is such that the final velocity is much 
smaller than the speed of sound c s . The aim was to check whether the system is Galilean 
invariant and the drop is correctly convected by the flow. We monitored the shape of the 
drop and measured its center of mass velocity v CM . This is defined as the average velocity 
of the center of mass whose position is 

TOM® = E y jMt) (40) 
11 



where the sum is over the lattice nodes r^- inside the drop. This velocity represents the 
convection velocity and is compared with the fluid velocity V/(i) = u*(r cM(t)) at the center 
of mass given directly by the LBM. In Fig. 5 the comparison between the velocities vcm 
and \f along the x-direction is shown in the case with g x — 3 x 10~ 5 . It is evident that the 
two coincide indicating that the drop is correctly advected by the fluid. Moreover, its shape 
is not altered by motion as it can be seen in Fig. 6 where some configurations of the system 
at different times are presented. Moreover, the drop is shown to make clear that it does not 
change in shape with time. We measured the ratio of the horizontal and vertical diameters 
finding that it stays almost constant with a deviation less than 3% from the value 1. If the 
advection velocity is higher, the drop will be slightly deformed being stretched along the 
x-direction. This effect becomes negligible when increasing the surface tension (J3J) via the 
parameter k. 

IV. CONCLUSIONS 

In this paper we have considered a lattice Boltzmann method for binary mixtures with 
thermodynamics fixed by a free-energy functional. We used a mixed method, with conti- 
nuity and Navier-Stokes equations simulated by LBM, and convection-diffusion equation by 
finite difference schemes. Differently than in previous free-energy LBM formulations [71], the 
interaction part in the pressure tensor is not introduced by fixing the second moment of the 
LBM populations but by introducing a forcing term in the lattice equation. This approach is 
suggested by a microscopic picture and allows to obtain a continuum limit without spurious 
terms. On the other hand, the mixed or hybrid approach allows a reduction in the required 
memory and this can be relevant in performing large-scale simulations. 

In order to reduce spurious velocities, differential operators have been discretized by gen- 
eralizing the usual lattice representations. Free parameters appear and their optimal values 
have been fixed by requiring that the maximum value of spurious velocities at equilibrium 
is minimized. 

We considered simple test situations, flat interfaces and single drops showing that the 
correct equilibrium profiles are reproduced. We found that spurious velocities are reduced 
of about an order of magnitude when a more general stencil is applied to the derivatives 
in the forcing term of the LBM equations. We did not found any relevant difference by 
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applying this procedure to the differential operators appearing in the convection-diffusion 
equation. We also checked that our method is stable in phase separation studies, even if 
we have not reported the results of these simulations in this work. Finally, we checked the 
effective Galilean invariance of the system by advecting for some time interval by a constant 
force a configuration with one drop and then letting the system to evolve without forcing. 
For the cases considered, we did not observe relevant drop deformations, the drop being 
correctly advected by the surrounding fluid. 

In conclusion, we hope that this development of the free-energy LBM can be useful in 
future simulations of binary mixtures and complex fluids. 
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FIG. 1: Cell of the D2Q9 lattice used in the present study. 
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FIG. 2: Equilibrium profile of a planar interface on a lattice of size L x = L v = 64 in the SC case. 



The continuous line is the analytical result (|2|) and data points are the results of simulations. 
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FIG. 3: Velocity patterns (the same scale is used in both the panels) at equilibrium when r/AtLB = 
5 in the SC case (upper panel) and in the OC case (lower panel). Empty spaces are due to negligible 
values of velocity. In both the cases the system has size Lx — Ly — 128. 



17 




18 



M 0.02 
o 

x 

H — 



w 0.015 - 
o 




1 I I I I I I I I I I I I I I L 



2000 4000 6000 8000 

t / At LB 

FIG. 5: Velocities of the center of mass of the drop vqmx (°) and of the fluid vf x ( ) at the 

center of mass along the x-direction as a function of time. The external force G x acts until the 
time t/At LB = 500. 
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FIG. 6: Configurations of the advected drop at consecutive times. The system has size L x = 
256, L y = 128. In the lower panel the drop, extracted from the system, is shown on an underlying 
mesh to better appreciate its shape. 



20 



